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Abstract This paper extends the Method of Partic- 
ular Solutions (MPS) to the computation of eigenfre- 
quencies and eigenmodes of plates. Specific approxima- 
tion schemes are developed, with plane waves (MPS- 
PW) or Fourier-Bessel functions (MPS-FB). This frame- 
work also requires a suitable formulation of the bound- 
ary conditions. Numerical tests, on two plates with var- 
ious boundary conditions, demonstrate that the pro- 
posed approach provides competitive results with stan- 
dard numerical schemes such as the Finite Element 
Method, at reduced complexity, and with large flexi- 
bility in the implementation choices. 

Keywords clement-free methods • biharmonic 
equation • numerical methods • algorithms • eigenvalues 

1 Introduction 

Numerical computation of eigenfrequencies and eigen- 
modes of plates is an important problem in mechanics. 
An eigenmode is a non-zero solution to 
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with boundary conditions, where D is the rigidity of the 
plate, p the specific mass of its material, h its thickness, 
and T the normal tension applied at its edges, assumed 
to be uniform. The eigenfrequencies are the frequencies 
Lo such that a non-zero solution exists. Apart from par- 
ticular cases where these quantities can be analytically 
computed (e.g., circular plates with simple boundary 
conditions), they must be obtained by numerical meth- 
ods. 

The finite element method (FEM) , which uses piece- 
wise polynomials to approximate the solutions, can be 
used to compute these eigenmodes and eigenfrequen- 
cies. However, it can be computationally intensive at 
high frequencies, as the size of the numerical problem 
scales as the square of the spatial frequency. Alterna- 
tive methods are the boundary element method (BEM) 
[1], the method of fundamental solution (MFS) [2, or 
its variant proposed by Kang et al. [3], the Non Di- 
mensional Influence Function (NDIF). Here, only the 
solutions to the equation with a given wave number 
are approximated as linear combinations of a family of 
functions. Eigenmodes are found as such combinations 
which also satisfy the boundary conditions. 

Here, we investigate a new computational method 
derived from the Method of Particular Solutions, pro- 
posed by Fox, Henrici and MoUer (FHM) [1], and im- 
proved by Betcke and Trefethen [5] , for the computation 
of eigenmodes of the Laplace operator. While the fun- 
damental ideas used in this method are similar to the 
previously cited methods, it has several advantages: 

— the stability of the numerical problems is improved, 

— multiple eigenvalues are easier to determine, 

— and basis of the associated eigenspaces are readily 
available. 
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The goal of this paper is to extend this MPS frame- 
work, presented in section 2, to the computation of 
eigenmodes and eigenvalues of plates. Because of fun- 
damental limitations of the Vekua theory, we apply the 
method only on star-shaped plates with smooth bound- 
aries. Our first contribution is the analysis of an ap- 
proximation scheme based on the Vekua theory, given 
in section 3. This provides some bounds on the approxi- 
mation error of a solution of eqn. ([l]) by sums of Fourier- 
Bessel functions in Sobolev norms, based on similar re- 
sults for the simpler case of the Helmholtz equation. 
Our second contribution, described in section 4, is the 
formulation of the problem in a way compatible with 
the MPS, and its numerical evaluation presented in sec- 
tion 5. We discuss the extension to more general shapes, 
and various implementation matters, in section 6. 



pute its determinant d{k): 



2 The method of particular solutions 

The method of particular solutions was introduced by 
Fox, Henrici and MoUer (FHM) [3] for the case of eigen- 
modes of the Laplace operator in a L-shaped domain 
with a singular corner. 

The basic idea of this method is, instead of con- 
sidering the entire space in which the eigenmodes are 
searched (e.g. for the case of the Laplace opera- 
tor, approximated by finite element spaces), to con- 
sider separately the spaces of solutions of the Helmholtz 
equation for different wave numbers. Then, in each of 
these spaces, we can look for a nonzero function which 
also satisfies the boundary conditions, i.e., an eigen- 
mode. While building approximations for a lot of dif- 
ferent spaces seem to be counter-productive compared 
to the unique approximation needed for a Galerkin ap- 
proximation, this alternative scheme is interesting as ef- 
ficient approximations can be obtained for these spaces 
of solutions to the Helmholtz equation. 

The method developed by FHM, for the computa- 
tion of eigenmodes of the Laplace operator with Dirich- 
let boundary conditions, is as follows. For each wavenum- 
ber fc, we consider N points Xj on the border of the 
domain, and a family of N functions (jji spanning a 
subspace which approximates the set of solutions to 
the Helmholtz equation. The considered family, Fourier- 
Bessel functions of fractional orders, was specifically 
constructed to take into account the singularity aris- 
ing in the reentrant corner of the L-shaped domain. In 
order to find the eigenfrequencies, one has to construct 
a square matrix M{k) that contains the values of the 
functions at the N points of the border, and to com- 



d{k) = detM(fc) 



h{xi) ■ ■ ■ (t>N{xi) 
blixN) ■ ■ ■ (j)N{xN) 



(2) 



If k is an eigenfrequency, there is a non-zero solution to 
the Helmholtz equation with values zero on the bound- 
ary, and its approximation ^ ai4>i is thus close to zero 
at the sampling points, with nonzero coefficients a^. 
The image of the vector of coefficients (ai) by the ma- 
trix M{k) is precisely the vector of the values of ^ ai4>i 
on the points of the border. This means that the deter- 
minant of the matrix M{k) is close to zero. Therefore, 
the eigenfrequencies are obtained as local minima of 
d{k). 

Numerous variants of this method, using different 
approximation schemes, have been developed since the 
original article. They mostly difi^er by the functions used 
to approximate the solutions: the MFS uses fundamen- 
tal solutions, the NDIF [3j uses Bessel functions of the 
first kind of order 0, etc. 

As pointed out in [5 , this simple method has known 
limitations. The discretization of the space of solutions 
and the sampling of the boundaries must be of the same 
size, and more importantly, the matrix M (k) gets ill- 
conditioned as the number N of functions grows. In- 
deed, using a larger family of functions (/>i, in order to 
have better approximations of the modes, makes the 
problem numerically unstable. An interpretation of this 
fact is that the algorithm does not search for a non- 
zero function inside the domain with value zero on its 
boundary, but actually for a function with non-zero co- 
efficients of its expansion, with zero value on the bound- 
ary. The properties of the approximating families are 
such that having non-zero expansion coefficients does 
not ensure significant values of the function inside the 
domain. 

To avoid these problems, Betcke and Trefethen [5] 
suggest to solve, for each fc, the following optimization 
problem : 



T(fc) = min ||ru||^ 



^ such that 



{an) 



ML2{n} = 1 i^nd Au+k 



(3) 



where T is the trace operator on the boundary dO. 
Then, k is an eigenfrequency if and only if t(A:), called 
the tension, is zero. 

This problem can be discretized as follows. A family 
of functions ((pi) is chosen for the approximation of the 
solutions of the Helmholtz equation. For a function with 
expansion coefficients u = {ui, . . . , mat). 



I ^^11 is (ar2) 



u*Fu 
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u*Gu 



where the coefRcients of the matrices F and G are 



(4) 



These scalar products can be estimated by samphng 
the domain and its border, and using a Monte-Carlo 
approximation of the integrals. The optimization prob- 
lem (H can be replaced by the generalized eigenvalue 
problem 



AFu ^ Gu 



(5) 



for which the largest eigenvalue is the inverse of T{k). 
Note that here, the size of the matrices is the size of the 
approximating family, and does not depend on the num- 
ber of samples used in the domain and on its boundary. 

The eigenfrequencies are found as the local minima 
of T(fc). This method offers additional advantages. The 
coefficients of the expansion of the eigenmodes are read- 
ily available as the first eigenvector of eqn. ^ , and n- 
multiple eigenfrequencies are characterized by the fact 
that the n first eigenvalues of eqn. ^ exhibit a local 
minima. Here, a basis of the eigenspace is obtained us- 
ing the n first eigenvectors of eqn. Note that this 
basis has no reason to be orthogonal, as it is the vectors 
of the coefficients of the expansions of the basis func- 
tions that are orthogonal, not the functions themselves. 

The same method, with a slightly different imple- 
mentation, was used by Barnctt and Berry [S] to com- 
pute high frequency modes of quantum cavities. As they 
considered only domains with smooth boundaries, plane 
wave families were sufficient for the application of the 
method. 

To generalize this method to plates, two adaptations 
are needed: 

— an approximation scheme for solutions of eqn. ([T]) 
has to be developed 

— the boundary conditions encountered in plate prob- 
lems have to be modeled in a way compatible with 
formulation eqn. ([S]). 

These two points are the topics of the next two sections, 
respectively. 



3 Approximation of plate eigenmodes 

In this section, we prove that solutions of eqn. ([T]), with 
arbitrary boundary conditions, can be approximated by 
sums of Fourier-Bessel functions and modified Fourier- 
Bessel functions. We first give a short account of the 
Vekua theory for the Laplace operator, and then extend 
these results to the bi-Laplace operator. 




Fig. 1 A star-shaped domain satisfying the exterior cone 
condition with angle Att 



The domain where the functions are defined is as- 
sumed to be star-shaped, and to contain the ball cen- 
tered on the origin of radius ph, where h is the diam- 
eter of the domain. Furthermore, we assume that the 
domain satisfies the exterior cone condition with angle 
Att. This means that each point of the border is the 
vertex of a cone of angle Att which does not intersect 
the interior of the domain. Such a domain is pictured 
on figure [TJ 

In the following, we give the approximation bounds 
in weighted Sobolev norms defined by 



\m,k 



P=0 Pl+P2=P 



QP 



dxdy (6) 



3.1 Vekua theory for the Laplace operator 

A simple example of an approximation of a solution 
to a differential equation is the case of holomorphic 
functions. These functions, solutions to the Cauchy- 
Riemann equations, can be, on a simply connected do- 
main, approximated by polynomials of the complex vari- 
able. By taking the real part of an holomorphic function 
and of its approximations, it is shown that we can ap- 
proximate an harmonic function, solution to Au — 0, 
by harmonic polynomials in R^. 

The Vekua theory, exposed in f?^ and summarized in 
[8] , gives similar results for solutions of elliptic PDEs by 
generalizing the operation "taking the real part", which 
allows us to map holomorphic functions to harmonic 
functions, to solutions of these PDEs. Using these oper- 
ators, which are continuous and continuously invertible, 
approximation of holomorphic functions by polynomials 
of the complex variable is translated to approximation 
of solutions to the PDEs by the images of polynomials. 
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In the case of the Helmhohz equation 



Furthermore, if w ^ H^^^, th 



en 



Au + \u = {) 



(7) 



Further details on these operators and their properties 
in Sobolev spaces can be found in [3]. 

The main result of Moiola et al. is that the solu- 
tions of eqn. ([?]) can be approximated by generalized 
harmonic polynomials, i.e. functions of the form 



Q;„J„(fcr)e* 

n— — L 



(8) 



where (r, 9) are the polar coordinates. 

When A is strictly positive, Moiola et al. [S] have 
shown that there exists a generalized harmonic polyno- 
mial of order at most L such that 



\\u - QlWj.m < Cil + (fc;i)^+6)e3(i-'')'=V4 



log(^ + 2) 
L + 2 



\{K+l-j) 



(9) 



m\K+i,k 



In the case where A < 0, the bound is multiplied by 

g3fch/2 

When the function to be approximated is infinitely 
differentiable in an open domain containing f?, the con- 
vergence is exponential in L jlOj . 

Moiola et al. used this results to analyze the ap- 
proximation of solutions of eqn. ([7| by sums of plane 
waves [11]. The orders of convergence are identical to 
the generalized harmonic polynomials case. 



3.2 Extension to plates eigenmodes 

Eigenmodes of an homogenous plate of rigidity D, spe- 
cific mass p and thickness h, with in-plane tension T, 
are solutions of 



DA^w + TAw - phuj^w 0. 



(10) 



In order to approximate such functions, we reduce 
this problem to the approximation of two solutions of 
the Helmholtz equation. Indeed, solutions of equation 



( 10 ) can be decomposed as a sum of two solutions of 



equation ([t]), with parameters deduced from the prop- 
erties of the plate. 



Lemma 1 Let w a solution of in the sense of dis- 
tributions. Then w can be decomposed as the sum of 
wi solution of Awi — XiWi = 0, and W2 solution of 
A'W2 — X2W2 — 0, where Ai and A2 are the zeros of 



DX^ + TX- phuj^. 



(11) 



2fc2 

\wi\\kM < -J^\\w\\k+2M 
\w2\\K,k2 < -^\\w\\k+2M2- 



(12) 
(13) 



where ki = y'\Xi\ and k2 = y'\X2\, being the largest, 
and Sx — \/T^ + 'iDphuj'^/D the difference between Xi 
and X2 

Proof Analysis : assuming the decomposition, we find 

Aw — X2W = {Awi — X2W1) + {Aw2 — X2W2) (14) 
= (Ai-A2)wi+0 (15) 

A similar computation for W2 gives 

1 , . , . 1 



Wi 



Ai-A: 



-{AW — X2W), W2 



X2 — Al 



(Aw — Xiw). 



Note that Ai and A2 are always distinct as Sx is always 
strictly positive. 

Synthesis : we check that w = wi + W2 ■ 



Wi + W2 



1 



Al — A2 
w 



(Aw - X2W) - {Aw - Xiw) (16) 

(17) 



Then that Awi — XiWi = : 
1 



Awi — XiWi — 



Al — A2 
1 

Al — A2 




{{A'^w - X2AW) - {XiAw - A1A2 

(18) 

(A'^w - (A2 + Xi)Aw + X1X2W) 

(19) 
(20) 



The last equality comes from the fact that Ai and A2 
are the zeros of the polynomial DX^ + TX — phui"^. We 
also find that Aw2 — A2U'2 = 0. 
Finally, if w g then 



ll'^ilk.fci < 



1 



< 



Al- A2I 
1 

Sx 
2kl 



{\\Aw\\K,k,+^2\\w\\KM) (21) 



< i: {kl\\w\\K+2M + ki\\w\\K+2M) (22) 
\Mk+2M (23) 



3A 



The result is identical for W2. 

Remark 2 In the case where no tension is applied to 
the plate, i.e. T — 0, we have ki — k2 — {ph/ DY^'^y^ = 
k, Sx — 2/c^ and the results can be simplified as : 

\\wi\\K,k < \\w\\K+2,k, ||lf2||i^,fc < 1 1 It; 1 1 if +2, fc- 
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Remark 3 Two orders are lost in the majorizations 
{ 12 ) and ( 13 ) ; the norm of order K of both compo- 
nents of w are hounded by their norm of order K + 2. 
It is impossible, in the general case, to have a better 
bound. Let consider, on a disk sector centered at the 
origin, the function f defined by f = e*^/^( Ji/2('') ~ 
Ii/2{'f)) in polar coordinates. This function is solution 
of A^w — w = 0, and can be decomposed as the sum 
of fi = e'^/^ Ji/2(f ), solution of Afi + fi = 0, and 
/2 = -e'''^^Ii/2{r), solution of Af2 - /2 = 0. These 
two functions have a radial behavior at the origin simi- 
lar to r^/^, and thus are not in . Their sum however 
behaves like r^/^, and is in H^. 

If Ai, or A2, is negative, the corresponding compo- 
nent of w can be readily approximated by generalized 
harmonic polynomials or plane waves using the results 
of Moiola et al. llj. If Ai or A2 is positive, the as- 
sociated component can be approximated in a similar 
way. In that case, we use either a family of modified 
Fourier-Bessel functions; where the Bessel functions J„ 
are replaced by modified Bessel functions /„, or a fam- 
ily of exponential functions e"^'^ instead of plane waves. 
The bounds on the approximation error are similar. 



Theorem 4 Let Q he a domain satisfying the assump- 
tions of this section, K > I integer, and w € ver- 
ifying conditions of lemma^ Then for all L > K, there 
exist two generalized harmonic polynomials and Ql 
with parameters Ai and A2, of degree at most L such 
that for all j < K , 



\w - [Pl 



Ox 



\n{L 



L + 2 



{k+h)''-^\w\\K+2,k^.. 

(24) 



where k- is the smallest of ki and k2, and fc+ the 
largest. If Ai and A2 are both negative, the bound can 
be divided by 



Proof The first step of the proof is to decompose w 
using lemmajl] Using theorem 2.2.1.ii and remark 1.2.6 
from |12j . we can approximate these two components 
by generalized harmonic polynomials, modified if A is 
positive. 

Let us assume Ai < 0, A2 > 0, with IA2I > |Ai|. 
Other cases can be treated similarly. 



We have 

\w - {Pl + (9l)||j,/c+ < llwillj^fc^ + \\w2\\j,k^ 

< \\wi\\]M + \\W2\\]M 

< C{l + {kihy+^)ei'^^-P^''''' 



ln(L + 2) 



HK-j) 



{k,hf-^\wi\\KM 



L + 2 

+ Cil + {k2hy+^)e^'^^-P^''''' 



HK-j) 



(k2h)''-^\wi\\K.M 



ln(L + 2) 
L + 2 

T 2 

< C-^(l + {k+hy+^)e-^^^-''^''+'' 



Remark 5 For plates without in-plane tension, the re- 
sult is slightly simpler. In that case, the roots of eqn. 



11) have same absolute values and opposite signs. We 



thus have 

\\w - [Pl + QL)Uk < C{1 + (A;/i)J"+6)ei(3-p)fc'' 
/ln(L + 2)^^(^"^^ 



{kh)''-'\\w\\K+2,k. 

(25) 



Remark 6 We prove the approximation result for sums 
of Fourier-Bessel functions. However, as Fourier-Bessel 
functions can he approximated by sums of plane waves, 
modified Fourier-Bessel functions can he approximated 
by sums of exponential function of the type e^ '^ with 
constant norm of k. This allows an approximation of 



solutions of eqn. ( 10 1 by sums of plane waves and ex- 
ponential functions. 



4 Boundary conditions 

Boundary conditions usually encountered in plate prob- 
lems (clamped edges, simply supported edges and free 
edges) are not as such readily usable for the numeri- 
cal scheme proposed here, and thus need to be modeled 
differently. 

For clamped edges, the displacement and its normal 
derivative are zero. The tension for clamped boundary 
conditions is then: 



tr. = 



1 

fc2 



dw 



dr 



A simply supported edge has zero displacement and tor- 
sion moment M„ 



t. = 



1 



D^k'i 



|M„| 
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Finally, free edges have zero torsion moment and Kelvin- 
Kirchhoff edge reaction Kn 



1 



\K„ 



The torsion moment and the Kelvin-Kirchhoff edge re- 
action writes [T3l 



-d(^ 



at2 



-D 



dn 



+ {2-1,) 



Q3 



W 



dndt^ 



B 



\ d'n? 



where u is the Poisson ratio of the material, and R the 
curvature radius of the boundary. 

The constants l/fc^, and in front of 

some integrals not only make the quantities homoge- 
neous, but also improve the numerical stability as the 
contribution of the two integrals are rescaled to have 
the same order of magnitude. For instance, in the case 
of tc estimated using plane waves, the first term con- 
tains products of planes waves, while the second term 
contains products of derivatives of plane waves, which 
are the plane waves themselves multiplied by the scalar 
product of the wave vector and a unit vector normal 
to the boundary. With no rescaling, the second term 
would be more and more influent as the frequency in- 
creases. This rescaling is similar to the weights used to 
define the norms || • \\m,k in eqn. (|6|. 

Note that for a plate with various boundary con- 
ditions along the border, the corresponding tension is 
the sum of the tension for the different boundary con- 
ditions, integrated on their respective domain. 

5 Numerical results 

We now compare the methods with analytical results for 
simple cases, and with numerical results obtained with 
Cast3M [14], a widely-used FEM simulation program. 
To avoid technicalities, the simulated plates are star- 
shaped with smooth boundaries. The treatment of more 
general shapes is discussed in the next section. Since 
these numerical tests are done without in-plane tension, 
we can use the wavenumber k — {p/ D)^/'^^/lo to express 
the eigenfrequencies. In the case of in-plane tension, the 
wavenumbers used to generate the Fourier-Bessel func- 
tions (resp. plane waves) and modified Fourier-Bessel 
functions (resp. exponential functions) are computed 
according to lemma [T] 

We first test the method on a circular plate of radius 
1, with various boundary conditions. Here, we use plane 
waves, as eigenmodes of circular plates are sums of a 
Fourier-Bessel function and a modified Fourier-Bessel 




Fig. 2 Inverses of the four largest eigenvalues (the first one 
being the tension) of discrete problem ijSJ, for the clamped 
circular plate 





Clamped 
Leissa MPS-PW 


Simply supported 
Leissa MPS-PW 


Leissa 


Free 
MPS-PW 


1 


3.196 


3.196 


2.23 


2.22 


2.29 


2.29 


2 


4.611 


4.611 


3.73 


3.73 


3.01 


3.01 


3 


5.906 


5.905 


5.06 


5.06 


3.50 


3.50 


4 


6.306 


6.306 


5.46 


5.45 


4.53 


4.53 


5 


7.144 


7.144 


n/a 


6.32 


4.65 


4.64 



Table 1 Low eigenfrequencies for a circular plate with var- 
ious boundary conditions, in Leissa |15| and computed with 
the MPS and plane waves (MPS-PW). 



function, making the numerical problem trivial. Table 
[T] gives the computed eigenfrequencies and the values 
given in Leissa [15^ for low-frequency modes. Figure [2] 
shows the tension as a function of the frequency in the 
clamped case, and the next three eigenvalues of the nu- 
merical problem ([s]) . These can be used to identify mul- 
tiple eigenvalues and to compute a basis of the associ- 
ated eigenspaces. Examples of eigenmodes are given for 
the three boundary conditions on figure [3] In these nu- 
merical experiments, the boundary of the disk was dis- 
cretized with 2048 points, its interior with 1024 points 
drawn using the uniform distribution on I?, and the 
number of plane waves (which is also the size of the 
numerical problems to solve) was taken as 8/c, which in 
this case ranges from 28 to 80. 

We now compute eigenfrequencies and eigenmodes 
of a second plate, of mode complex shape, with bound- 
aries defined by the parametric equations 



X = cos t 



2t 



y — s\Ta.t + sin ^ 



t e [0, 27r 



(26) 
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T3 
CU 
Q. 

E 
u 




# 




k=8.347 



k=3.196 



■o 

<D 

■M 
i— 

o 

Q. 
Q. 

in 

£ 
1/1 




k=5.451 



CU 




k=6.206 

Fig. 3 Some examples of eigenmodes, for simple and double eigenvalues, of the circular plate with clamped, simply supported 
and free boundary conditions 



Table [2] gives the first ten eigenfrequencies of this 
plate vifith the three types of boundary conditions : 
clamped, simply supported and free. The proposed method, 
MPS with plane waves (MPS-FB), uses 2048 points 
on the border and 1024 points inside. The number of 
planes waves is lOfc, which here ranges from 20 to 80. 
For the clamped conditions, the results using Fourier- 
Bessel functions are also given (method MPS-FB). We 
compare our results with those obtained by the FEM, 
as implemented in the CastSM package, with 6624 ele- 
ments and the border discretized by 180 segments. 

For two eigenmodes of the clamped plate, we give in 
table [3] the estimated eigenfrequencies with varying size 



of the discrete problems for FEM and the MPS-PW, 
and compare them to the values obtained in [2] with 
the MFS (with an unspecified size). These two modes, 
computed with the MPS-PW, are shown on figure |4] 



6 Discussion 

We now discuss some implementation details, relative 
to the treatment of more general shapes, or the accel- 
eration of the computations. 
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Clamped 




Simply supp. 




Free 




FEM 


MPS-PW 


MPS-FB 


FEM 


MPS-PW 


FEM 


MPS-PW 


1 


3.3498 


3.3782 


3.3562 


2.3819 


2.4061 


2.2189 


2.2000 


2 


4.5931 


4.5944 


4.5864 


3.6789 


3.6733 


2.3473 


2.4001 


3 


4.8717 


4.9275 


4.8815 


3.9539 


3.9993 


2.7643 


2.7481 


4 


5.8588 


5.8376 


5.8366 


4.9799 


4.9505 


3.4544 


3.4222 


5 


6.1655 


6.2067 


6.1987 


5.2874 


5.3216 


3.7476 


3.7623 


6 


6.3869 


6.4567 


6.4527 


5.4832 


5.5456 


4.1207 


4.0993 


7 


7.1179 


7.0738 


7.0738 


6.2658 


6.2147 


4.2170 


4.2004 


8 


7.4596 


7.4949 


7.4949 


6.5912 


6.6188 


4.8551 


4.7985 


9 


7.7190 


7.7740 


7.7590 


6.8597 


6.9058 


4.9701 


4.9845 


10 


7.8910 


7.9300 


7.9660 


6.9979 


7.0308 


5.4432 


5.4316 



Table 2 First ten eigenfrequencies, expressed in wavenumber k, of the second plate with various boundary conditions, for 
FEM, MPS with plane waves (MPS-PW) and MPS with Fourier-Bessel functions (MPS-FB, clamped conditions only) 



MFS 


FEM 


MPS-PW 




246 


number of elements 
1018 6624 


26966 


number of plane 
50 60 


waves 
70 


9.11259 
9.27903 


8.9329 
9.1025 


9.0245 9.0590 
9.1830 9.2122 


9.0637 
9.2165 


9.1068 9.1121 
9.2651 9.2722 


9.1126 
9.2787 



Table 3 Eigenfrequencies for two modes of the clamped second plate, with various discretization sizes, computed by FEM (as 
implemented in Cast3M) and MPS-PW. Comparison with MFS (2j 




Fig. 4 Eigenmodes of the second plate, for k = 9.1126 
(left) and k = 9.2787 (right) computed by MFS-PW, for the 
clamped boundary conditions 



6.1 Fourier-Bessel vs. plane waves 

As shown by the numerical experiments, both Fourier- 
Bessel functions and plane waves can be used to approx- 
imate eigenmodes. They have similar approximation 
properties, but differ implementation-wise. Fourier-Bessel 
functions are orthogonal on a disc, ensuring better sta- 
bility, while plane waves are more and more ill-conditioned 
as their number increases. However, this can be treated 
by pre-conditioning the plane waves family with a dis- 
crete Fourier Transform, mapping the plane waves to 
approximations of the Fourier-Bessel functions. 

The main advantage of the plane waves is the straight- 
forward computation of their derivatives : differentiat- 
ing a plane wave along a certain direction amounts to 
multiplying it with the scalar product of its wave vector 
with a unit vector. This makes the construction of the 
matrices both easy to implement and fast. 



6.2 Shapes 



The proposed method relies on an approximation scheme 
for solutions to the equation ([l]), and therefore is lim- 
ited to cases where such approximations are available. 
In particular, the approximation scheme assumes a star- 
shaped domain. In the case of a simply connected, but 
non-star convex, domain, the approximation is not guar- 
anteed to succeed. A possible way to overcome the prob- 
lem is to cut the domain into star-convex subdomains, 
to approximate the solutions of ([T]) in these subdo- 
main, and to add terms in the tension ensuring that the 
displacement, normal derivative of the displacement, 
bending moment and strain have the same value at both 
sides of the internal boundaries. A similar method as al- 
ready been applied to the particular case of polygonal 
membranes 1161 ■ 



Other cases of non-convex domains are domains with 
holes, but star-convex if the holes are filled. In that 
case, following Vekua, one can approximate solutions of 
([!]) by adding to the family of planes waves or Fourier- 
Bessel functions, the sets of Fourier-Bessel functions of 
second kind e™^r„(fcr) and e*"^ii:„(fcr), one for each 
hole, centered on a point chosen in each of them. This 
type of approximation can be compared to the approx- 
imation of holomorphic functions given by the Runge 
theorem. An application to membranes can be found in 



Title Suppressed Due to Excessive Length 
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6.3 Singularities 

The domains considered here have smooth boundaries. 
This guarantees that the eigenmodes of the plates are 
smooth, and that the convergence of the approxima- 
tions is fast. However, as shown in theorem |4] singulari- 
ties, which can appears at corners of a domain with non- 
smooth boundaries, slow down the convergence of the 
approximations. The size of the discretization for such 
cases is then larger than for smooth boundaries, possi- 
bly too large to guarantee the numerical stability of the 
computations. To accelerate the convergence, FHM and 
BT use fractional Fourier-Bessel functions centered on 
the singular corner of the domain (in their numerical 
experiments, an L-shaped polygon). This idea as been 
used for plates in a different setting by De Smet et al. 

Note that, for a plate with polygonal holes, at least 
a corner of the hole is singular. In that case, it is im- 
possible to use fractional Fourier-Bessel functions, as it 
is impossible to define them on a domain containing a 
path around the origin. In that case, one can combine 
fractional Fourier-Bessel function with the method de- 
scribed in the previous section [T7] . 



6.4 Numerical considerations 

Although the numerical stability is improved compared 
to the determinant based MPS, the improved version 
is still prone to instabilities when a large approxima- 
tion order is used. This instabilities are further am- 
plified in the case of plates, where modified Fourier- 
Bessel functions, or exponential functions, are included 
in the approximating families. The behavior of these 
functions are such that they are non-negligible only on 
a small region near the boundary of the domain. Using 
the Monte-Carlo approximation with uniform density 
to estimate the coefficients of the matrices is thus un- 
stable. Using a non-uniform density of samples, with 
more samples near the border would be a way to im- 
prove the stability of the estimations, in a way similar 
to what is used in [15] , where the reconstruction of a 
solution of the Helmholtz equation on a disc is improved 
by placing a fraction of the samples on the border of the 
disc. 



6.5 Speeding up the eigenvalue search 

In order to locate the minima of the tension, the pro- 
posed algorithm simply computed it on linearly spaced 
values in the interval we were interested in. However, 



the particular behavior of the tension (and more gener- 
ally, of the eigenvalues of problem ([s])) could be used to 
accelerate the search of these minima. Indeed, the ten- 
sion is a sequence of branches which behave more or less 
as parabolas. Newton iterations, along with the compu- 
tation of the derivatives of the eigenvalues of problem 
([5]), could therefore be used to quickly locate the min- 
ima. 



7 Conclusion 

This paper has described the extension of the Method 
of Particular Solutions for the computation of eigen- 
modes of plates. This method has numerous advan- 
tages. It can be used with any approximation scheme 
for the solutions of the studied equation ; in this paper, 
we used Fourier-Bessel functions and plane waves, but 
the method could be extended with fractional Fourier- 
Bessel functions in order to treat singularities. Its for- 
mulation also offers a large flexibility in the discretiza- 
tion of the domain, and independently, in the size of 
the numerical problem. The determination of multiple 
eigenvalues and eigenmodes is also straightforward. Fi- 
nally, the so-called tension, that has to be minimized to 
find the eigenfrequencies, has a specific shape that can 
be used to speed up the search. Future improvements 
include sampling schemes yielding better stability of the 
numerical problems, and using this method to compute 
high frequency eigenmodes of plates in an efficient way. 



Reproducible research 

The Matlab/Octave code used to compute eigenfrequen- 
cies and eigenmodes of figures [3] and [4] is available online 

m- 
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